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1.  INTRODUCTION 


Modeling  composite  materials  poses  a  significant  computational  challenge,  requiring  the  capability  to 
implement  the  material  response  accurately  and  efficiently  optimizing  the  computational  cost.  To 
accurately  simulate  composite  material,  the  orthotropic  effects  should  be  incorporated  by  importing 
empirical  data  [1]  or  using  assumptions  within  theoretical  calculations.  The  most  notable  contribution  to 
the  composite  material  properties  is  the  reaction  from  the  local  micro-structure  especially  in  critical 
regions,  such  as  high  temperature  changes.  However,  within  these  critical  regions  incorporating  the 
microstructure  effect  has  proven  to  be  problematic.  A  full  finite  element  solution  traversed  down  to  the 
microscopic  level,  taking  into  account  the  detailed  micro-structural  features  within  these  areas,  is 
currently  computationally  infeasible.  Hence  a  robust  and  reliable  multiscale  computational  technique  is 
necessary  to  account  the  microstructure  phenomena. 

Though  empirical  data  can  be  used  to  approximate  the  material  properties  for  an  ideal  non-damaged 
specimen,  performing  experimentation  on  damaged  composites  and  extracting  a  correlation  between 
micro-structure  effects  to  material  properties  can  be  extremely  difficult.  As  an  alternative,  material 
properties  from  a  composite  material  can  be  extracted  using  a  numerical  approximation  method,  referred 
to  as  the  homogenization  method  [2],  This  method  uses  asymptotic  expansions  of  field  variables  about 
macroscopic  values  and  provides  overall  effective  properties  as  well  as  microscopic  stress  and  strain 
values  for  a  local  point  on  the  global  model.  The  primary  limitation  with  this  method  is  that  it  assumes 
uniformity  of  the  macroscopic  fields  within  each  representative  volume  element  (RVE).  Hence,  this 
method  breaks  down  in  critical  regions  of  high  gradients.  To  resolve  the  limitation  of  the  homogenization 
method  within  critical  regions,  several  multiscale  methods  have  been  presented  which  are  reviewed  by 
Macri  et  al.  in  reference  [3], 

To  overcome  the  drawbacks  of  the  existing  methods  [3,  4]  proposed  a  structural  based  enrichment 
method  based  on  the  principles  of  partition  of  unity  which  allows  enrichment  of  the  approximation  space 
in  localized  sub  domains  using  specialized  functions  that  may  be  generated  based  on  a  priori  information 
regarding  asymptotic  expansions  of  local  stress  fields  and  microstructure.  One  of  the  major  advantages  of 
the  partition  of  unity-based  enrichment  strategies  is  that  the  enriched  local  function  space  may  be  easily 
varied  from  one  node  to  the  other. 

In  this  presentation,  we  focus  on  the  incorporation  of  non-linear  material  properties  into  the  simulation. 
Most  composites  will  be  composed  of  materials  that  are  temperature  dependent.  As  the  composite  is 
heated,  the  simulation  will  constantly  need  to  check  and  update  the  material  properties.  As  documented 
in  [5],  this  translates  to  updating  the  effective  properties  for  the  homogenization  method  at  every  time 
step.  In  [6],  a  similar  approach  is  addressed  for  structural  enrichment  occurring  in  mechanical  simulations. 
In  this  paper,  we  discuss  how  temperature  dependencies  are  implemented  for  thermal  simulations. 

In  the  next  section,  we  briefly  review  the  homogenization  and  enrichment  methods  and  how  temperature 
dependent  material  properties  are  incorporated  into  the  simulation.  And  in  section  3,  we  demonstrate  an 
example  of  the  approach. 
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2.  EXPERIMENTATION 


We  have  divided  this  section  into  two  parts.  Section  2.1  reviews  the  concepts  of  the  homogenization 
method  and  how  non-linearity  is  implemented.  In  section  2.2  we  discuss  how  this  is  implemented  for 
structural  enrichment. 


Figure  1.  Multi-scale  homogenization  approach 


2.1  Review  of  the  Homogenization  Approach 

The  homogenization  theory  [2],  assumes  multiple  scales  related  by  a  scale  factor,  y.  In  the  case  of  two 
scales: 


y 


X 

y 


[i] 


where  x  and  y  represent  a  global  and  microscale,  respectively  (see  Figure  1).  A  function  spanning  the 
scales  can  be  defined  by 

fY(x)  =  f{x,y(x))  [2] 

Such  that  the  microscale  is  dependent  on  the  global  scale.  The  homogenization  approach  assumes  local- 
periodicity  on  the  micro-scale.  Thus,  a  representative  volume  element  (RVE),  as  shown  in  Figure  1,  with  a 
side  length  of  Y can  be  depicted  as: 

f(xi,yj)  =  f(xi,yj  +  kYj)  [3] 

where  k  is  the  periodic  interval.  Using  asymptotic  expansion,  we  can  write  the  field  variables,  such  as 
displacement,  as: 


Ty(x,y )  =  +yT1(x,y)  +  0(y2) 


[4] 
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In  this  presentation  we  assume  a  thermo-elastic  response  for  the  system.  The  equations  governing  this 
thermal  analysis  are: 

)  +Q-  pcTY  =0  e  n  [5a] 

,Xi 

TY  =  Tg  e  rg  [5b] 

~kijTJcjni  =Cl  e  rq  [5c3 

The  k  term  represents  thermal  conductivity,  Q  is  a  heat  supply  per  unit  volume,  and  Tg  and  q  are  the 
prescribed  temperatures  and  heat  flux,  respectively. 

The  thermal  strains  are  calculated  such  that: 


£kl  —  aki9Y  ^ 

Where  dv  is  the  change  in  temperature,  defined  as: 

6y  =  ty  -  t,  [7] 


where  T,  is  the  initial  temperature.  The  change  in  temperature  field,  0Y,  will  be  used  to  calculate  the 
thermal  strains.  In  equation  [6],  a  is  the  coefficient  of  thermal  expansion  tensor. 


To  determine  the  effects  of  the  micro-structure  across  the  domain,  we  substitute  the  asymptotic 
expansion  of  the  temperature  into  the  heat  equation,  [5],  Rearranging  the  terms  gives  us: 


2  [O^7^),  J  + 

1  [(kiJT,xj)yi  +  {ktJT.yj)iXi  +  (/ci/7J'j)yJ  + 

0  [mo* + (ki^,yi + + q  ■  pcTi 


[8] 


+  0(y)  =  0 


By  setting  each  order  yl  term  to  zero,  it  is  straightforward  to  derive  the  following 
T\x,y)  =  T\x) 


[9a] 


[k.j  (lmj  +  <y;)]  ^  =  0  enRVE  [9b] 

[kijTx]  +Q~  pcT°t  =  0  en  [9c] 

,Xi 

Equation  [9a]  states  that  T°  is  only  dependent  on  the  global  scale,  x.  Equation  [9b]  is  the  governing 
equation  for  the  microscale,  where  /  is  the  identity  tensor  and  Hk(y)  is  a  Y-periodic  tensor  on  the  micro¬ 
scale  relating  T1  to  the  macro-scale,  such  that: 

T\x,y)  =  H^(y)T0Xm(x)  [12] 


where  m  is  based  on  the  dimensionality  of  the  problem.  Equation  [9c]  is  the  governing  equation  for  the 
global  scale.  In  this  equation,  k  represents  the  effective  thermal  conductivity  tensor,  p  is  the  effective 
density,  c  is  the  effective  specific  heat  and  Q  is  the  effective  heat  supply.  They  are  defined  as: 
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[13] 


For  an  implicit  time  simulations  the  well  documented  weak  form  of  9c  can  be  written  as 
(M  +  A tK)t+Ata  =  AtF  +  M  la 


[14] 

[15] 

[16] 


[17] 


-l 


M  =  N1  pcNdtt 


K  =  i 

J a 


BTkBd(l 


[18] 

[19] 


F=  [  NTQ(t  +  At)dn+  [  NTq(t  +  At)dT  [2°] 

■'n  Jr 

where  N  is  the  finite  element  shape  functions,  B  are  the  derivatives  of  N  with  respect  to  coordinates,  At 
is  the  change  in  time  and  q is  the  heat  flux.  The  global  temperature  at  time,  t,  is  defined  as: 

tT°(x )  =  N(x )  ta  [21] 


For  simplicity,  we  assume  problems  consisting  of  temperature  dependent  thermal  conductivity.  Thus, 
equation  [13]  becomes: 


^  H‘+AlT%  (/„;  +  a/J"™  1221 

For  linear  temperature  independent  simulations,  k  is  calculated  once  at  the  beginning  of  the  simulation. 
However,  since  it  is  assumed  that  t+AtTY  will  vary  across  the  model  and  with  each  time  step,  unlike  linear 
simulations,  k  will  need  to  be  calculated  at  each  integration  point  within  the  model,  for  each  time  step. 

For  an  implicit  simulation  we  calculate  k  at  time  step  t+At.  Thus,  equations  [17]  and  [19],  can  be  rewritten, 
respectfully,  as: 

(M  +  A tt+MK)  t+Ata  =  AtF  +  M  za  [23] 


t+MK  =  [  BTt+AtkBdn 


[24] 
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Because  of  K's  dependency  of  time,  equation  [23]  requires  a  non-linear  approach  at  each  time  step  to 
solve  for  t+Ata.  Using  Tayler  series  expansion  of  [23]  about  t+Ata  and  the  Newton-Raphson  method  we  can 
write  the  iterative  equations  as: 

j(j-VAaj  =  A tF  +  M  ta-  (M  +  Att+Atffy_1))t+Ataw_1)  [25] 


t+Ata(;')  =  Aa;  +  t+Atao-i) 


[26] 


where  /M'  is  the  derivative  of  [23]  with  respect  to  t+Ata  that  is  dependent  on  the  relationship  between  the 
thermal  conductivity  and  temperature. 

Thus  the  process  at  each  time  step  is  outlined  as: 

1.  Set  t+Ata°  =  *01  for  all  nodes  not  on  an  essential  boundary 

2.  Set  j=l 

3.  While  (not-converged) 

a.  Calculate  k(t+AtT^~1^  for  each  integration  point  within  the  model 

b.  Setup  Ju~^Aaj  =  A  tF  +  M  ta-  (M  +  A 

c.  Solve  for  Aa.1 

d.  Set t+Ata(j)=  t+Ata(i"1)  +  Aai 

e.  Set  j=j+l 

f ■  '/M,  <  tolerence,  set  status  to  converged 

4.  End  while  loop 

For  the  example  in  section  3,  we  will  be  examining  the  thermal  strains  generated  from  the  simulation.  In 
reference  [7]  we  define  the  effective  coefficient  of  thermal  expansion  as: 

&kl  =  yRVE  J  RVE  Cijkl(H<k,yl  -  ttkl)dflRVE 


2.2  Implementation  of  Damaged  Microstructure 

Within  this  section  we  review  the  partition  of  unity  paradigm  and  how  it  is  implemented,  followed  by  the 
derivation  of  the  enrichment  functions. 

2.2.1  Interpolation  Functions 

The  partition  of  unity  paradigm  requires  that  the  sum  of  all  functions  at  a  given  point  in  a  domain  is  equal 
to  one,  such  that 

N 

^?W  =  1  [28] 

i=i 

In  the  above  equation,  is  the  Ith  partition  of  unity  function  and  N  is  the  total  number  of  partition  of 
unity  functions  within  the  domain,  O.  Thus,  the  sum  of  all  the  functions  at  a  given  point,  x,  within  D  will 
equal  1. 
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Based  on  this  concept  we  define  the  global  approximation  space,  for  a  method  using  interpolation 
functions,  as: 

Vh,p  =  yN  (pOVh-P  c  //i(n)  [29] 

where  the  superscripts  h  and  p  are  the  size  of  the  element  and  the  polynomial  order,  respectively.  H1  is 
the  first  order  Hilbert  space  and  V [,l,p  is  the  local  approximation  space  at  /,  defined  as: 

V,h,p  =  spanme<{pm(x )}  c  H1(BIC 1X1)  [3°] 

In  equation  [17],  £  is  an  index  set,  B,  is  the  support  for  node  /  and  pm(x)  is  composed  of  polynomials  or 
other  functions.  From  the  above  equations  we  can  write  any  function  vh-p  within  the  global  approximation 
space  as: 

N 

vh-P(x)  =  ^  ^  hIm(x)aIm ,  vh’P  e  Vh’P  [31] 

1=1  me( 


where 


ft/mW  =  <P?(.x)Pm(.x)  [32] 

is  the  interpolation  function  at  node  /,  corresponding  to  the  mth  degree  of  freedom,  and  aIm  is  the 
associated  degree  of  freedom. 

For  FEA,  the  standard  shape  function,  /V,,  which  can  be  referenced  in  numerous  sources,  satisfies  the 
partition  of  unity  criteria  for  cp®.  For  a  typical  finite  element  analysis,  correlating  the  definition  of  a  finite 
element  shape  function  to  the  interpolation  function  in  equation  [32],  requires  that  pm  =  span{  1). 
Reference  [7]  discusses  other  partition  of  unity  functions  and  basis  that  can  be  implemented. 

2.2.2  Enrichment  Functions 

In  the  previous  section,  we  defined  a  function  in  the  global  approximation  space  as  the  sum  of  the 
partition  of  unity  functions  multiplied  by  a  group  of  basis  functions  and  associated  degrees  of  freedom. 
Focusing  on  the  basis  functions,  pn(x),  we  state  that  if  an  arbitrary  function  is  included  in  the  local  basis  of 
all  nodes  within  the  domain  and  assuming  aIm  =  Smn  V/  then 

•  N 

I 

1=1  met,  \/= 1 

Equation  [33]  demonstrates  that  it  is  possible  to  exactly  reproduce  a  function  over  the  entire  domain.  This 
fundamental  concept  enables  the  consistency  of  the  interpolation  functions  to  be  adjusted  based  on  the 
population  of  the  basis  functions.  Thus,  we  can  expend  the  basis  function  to  include  relevant  functions  to 
enforce  a  known  displacement  field. 

To  derive  the  structural  based  enrichment  function,  we  examine  the  asymptotic  expansion  of  the 
temperature  field.  Substituting  equation  [12]  into  [4]  gives  us: 

TY(x,y)  =  T\x)  +  y//£(y)r°n(*)  [34] 


II 


<P?(x)pm(x)aIrn  =  pn(x) 


0/W  =PnO) 


[33] 
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For  a  single-scale  problem,  we  define  the  displacement  field,  assuming  linear  consistency  in  2D,  as 

N 

T°(x)  =  11  0/°(*)Pm(*Mm  [35^ 

/= 1  me( 

where  and  pm  are  A//  and  pm  =  span{  1). 

To  perform  a  multi-scale  analysis,  we  need  to  incorporate  the  second  term  on  the  right  hand  side  of 
equations  [34],  We  focus  on  the  Y-periodic  tensors.  Details  of  the  derivation  are  given  in  [6],  The  resulting 
basis  functions  for  a  two-dimensional  thermal  problem  are  given  as: 

pm  =  span{  1,  Wf(y(x)),  H^(y(x))}  [36] 

Thus  using  equation  [32],  we  can  define  TY  as: 

TY(x)  =  hlm(x)alm  [37] 

where  <p° is  defined  as  /V,  and  pm  is  defined  by  equation  [36], 

2.2.3  Implementation  the  Multiscale  Enrichment  Technique 

The  approach  to  implementing  structural  based  enrichment  varies  depending  on  the  governing  method. 
In  this  presentation  we  will  focus  in  the  FEA  approach.  Reference  [4]  gives  complete  details  on  the 
implementation  of  the  multiscale  enrichment  into  FEA.  In  this  section  we  will  briefly  review  how  the 
algorithm  is  applied. 

The  analysis  is  based  on  a  homogenization-like  integration  (HU)  scheme  [4],  The  solution  requires  the 
enrichment  functions  and  the  material  properties,  which  are  only  known  on  the  micro-scale,  y.  To 
implement  this  technique  and  resolve  the  aforementioned  issue,  integration  will  be  performed  on  an  RVE 
transformed  to  the  macro-scale,  x,  coordinate  system.  The  HU  scheme,  shown  in  Figure  2,  is  based  on  the 
assumption  that  there  is  a  significant  deference  between  the  two  scales  and  that  the  Gauss  points  used 
for  FEA  are  sparse  with  respect  to  the  size  of  the  RVE.  In  this  approach,  an  RVE  is  centered  on  each  Gauss 
point  within  macro-scale  problem.  An  integrand,  /,  on  the  macro-scale  is  replaced  with  the  following 

nguass  Kguass 

I  =  j  WdSl=  £  WJiv(x?p)=  £  Wdi^WE  a^>RVE  [38] 

n  i= 1  i= 1  ni 

Where  -tr  is  an  arbitrary  function,  xEP  is  the  ith  Gauss  point,  HPVE  is  the  domain  of  an  RVE  centered  on  the 
ith  Gauss  point,  J  is  the  Jacobian  and  1/14  is  a  weight  associated  with  the  ith  Gauss  point.  Reference  [4]  has 
demonstrated  that  the  accuracy  of  HU  improves  with  decreasing  size  of  the  RVE.  Thus,  this 
implementation  works  well  when  the  difference  in  scales  is  significant.  However,  as  the  scales  become 
closer  and  the  therefore  the  volume  of  the  RVE  increases,  equation  [34]  will  develop  convergence  errors. 
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t'U)  tM,RVK 


Figure  2.  Implementation  of  multiscale  enrichment  into  FEA 


2.2.4  Implementation  of  the  Multiscale  Enrichment  Technique  for  Non-Linear  Simulation 

For  an  implicit  time  simulations  using  structural  enrichment  techniques,  the  weak  form  of  the  governing 
equation  is  setup  similar  to  the  homogenization  approach 


(M  +  AtK)t+Ata  =  A tF  +  M  la 

M  =  f  hT pchdfi. 

Jo. 

WkBdti 


K  =  l 

F=  I  hTQ(t  +  At)dfl  + 
■'n 


I  hT q{t 


+  At)ar 


[39] 

[40] 

[41] 

[42] 


where  h  is  defined  in  [37]  and  B  is  the  derivatives  of  h  with  respect  to  coordinates.  The  global  temperature 
at  time,  t,  is  defined  as: 


tT°(x )  =  h(x) 


[43] 


For  problems  consisting  of  temperature  dependent  thermal  conductivity,  equations  [39]  and  [41]  can  be 
rewritten,  respectfully,  as: 


(M  +  A  tt+AtK)t+Ata  =  A  tF  +  M 


t+AtK=  f 
Jo 


BTt+AtkBdQ 


[44] 

[45] 


Because  the  basis  of  h  is  composed  of  H we  need  to  examine  equation  [9b],  which  will  now  be  adjusted 
to  account  for  the  temperature  dependency. 


[k(t+AtTr)ij(lmj  +  H^yj)]y  =0  e/2 


RVE 


[46] 
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Since  it  is  assumed  that  t+AtTY  will  vary  across  the  model  and  with  each  time  step,  H £  will  need  to  be 
calculated  at  each  integration  point  within  the  model,  for  each  time  step. 

Using  Tayler  series  expansion  of  [44]  about  t+Ata  and  the  Newton-Raphson  method  we  can  write  the 
iterative  equations  as: 

Ju~1}Aaj  =  A tF^-V  +  Mu~^  za  -  (My_1)  +  [47] 

Where  J11'11  is  the  derivative  of  [44]  with  respect  to  t+ata  that  is  dependent  on  the  relationship  between  the 
thermal  conductivity  and  temperature. 

Thus  the  process  at  each  time  step  is  outlined  as: 

1.  Set  t+Ata°  =  *01  for  all  nodes  not  on  an  essential  boundary 

2.  Set  j=l 

3.  While  (not-converged) 

a.  Calculate  H ^  for  each  integration  point  within  the  model 

b.  Setup  eg  u  at  ion  [47] 

c.  Solve  for  AaJ 

d.  Set  t+Ata(j)=  t+Ata(H)  +  Aaj 

e.  Set  j=j+l 

/•  </IMI,  <  tolerence,  set  status  to  converged 

4.  End  while  loop 


3.  RESULTS 

We  examined  a  two-dimensional  cross  section  of  a  composite  wafer  with  a  hole,  depicted  in  Figure  3.  The 
dynamic  simulation  was  modeled  assuming  an  implicit  time  step.  A  temperature  change  on  the  hole  was 
implemented  to  induce  thermal  stress.  In  this  example,  we  examine  the  upper  right  quadrant  of  the  cross 
sectional  area.  The  composite  is  composed  of  AI2O3  fibers,  all  orientated  perpendicular  to  the  cross 
section,  within  an  Aluminum  matrix.  The  material  properties  are  given  in  Tables  1  through  3.  The  RVEs, 
were  modeled  using  150  quadratic  1st  order  elements.  The  boundary  conditions  on  the  global  model  were 
a  heat  source  that  linearly  increased  on  the  surface  of  the  hole  from  22°C  to  400°C  over  a  second.  The 
course  model  was  simulated  with  102  linear  triangular  elements.  The  results  were  compared  to  a  fine 
mesh,  Figure  3a,  composed  of  5.2x10s  quadratic  1st  order  elements  modeled  down  to  the  micro-scale. 

The  results  for  the  temperature  fields  across  the  model  are  shown  in  Figures  4.  The  results  show  a  similar 
temperature  field  after  Is  for  all  three  cases.  The  focus  of  our  examination  is  the  thermal  strain,  et, 
occurring  on  the  micro-structure.  Figures  5-7  show  the  thermal  strains  on  the  RVEs  at  five  locations  within 
the  model  and  visually  show  similar  results  for  all  three  cases.  However,  if  we  examine  the  error  in  the 
solutions  in  tables  4  and  5,  we  see  that  structural  enrichment  is  improving  the  simulation,  ranging  from 
an  error  reduction  of  0.5%  to  12%. 
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Table  1.  Material  properties  of  RVE 


Material 

a 

P 

c 

K1) 

(tonne/mm3) 

(mJ/tonne  sc) 

Aluminum 

2.2xl0'5 

2.7x10  s 

8.98x10s 

Al203 

4.9  xlO'5 

1.15xl0'9 

1.7xl09 

Table  2.  Thermal  conductivity  for  aluminum 


ec 

k  (mW/(mm?C)) 

22 

237 

77 

240 

126 

236 

227 

236 

327 

231 

527 

218 

Table  3.  Thermal  conductivity  for  Al203 


-C 

k  (mW/(mm?C)) 

92 

0.12 

204 

0.13 

315 

0.14 

404 

0.15 
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a)  Refined  Mesh 


b)  Course  Mesh 


Figure  3.  Composite  sample 
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a)  Refined  Mesh 


b)  Pure  Homogenization  c)  Pine  Structural  Enrichment 
Figure  4:  Temperature  fields  after  Is 
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Figure  6:  Thermal  strains  for  homogenization  after  Is 


13 


DISTRIBUTION  STATEMENT  A.  Approved  for  public  release;  distribution  is  unlimited. 


Figure  7:  Thermal  strains  for  structural  enrichment  after  Is 


Table  4:  Error  in  thermal  strains  at  gauss  points  of  course  mesh  for  homogenization 


Element  # 

Matrix  (Al) 

Fiber  (AI203) 

Avg  et  (mm/mm) 

Relative  Error  (%) 

Avg  et  (mm/mm) 

Relative  Error  (%) 

4 

1.381xl0'2 

3.7 

6.198xl0"3 

5.09 

45 

1.325xl0'2 

3.25 

5.951xl0"3 

4.81 

63 

5.454xl0'3 

3.97 

2.449xl0"3 

6.49 

99 

7.128xl0'3 

4.01 

3.200xl0"3 

5.92 

47 

1.327xl0'3 

10.01 

5.957xl0"4 

13.5 
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Table  5:  Error  in  thermal  strains  at  gauss  points  of  course  mesh  for  structural  enrichment 


Element  # 

Matrix  (Al) 

Fiber  (AI203) 

Avg  et  (mm/mm) 

Relative  Error  (%) 

Avg  et  (mm/mm) 

Relative  Error  (%) 

4 

1.377xl0'2 

3.48 

6.184xl03 

4.87 

45 

1.320xl0'2 

2.87 

5.927xl03 

4.43 

63 

5.264xl0'3 

0.5 

2.363xl03 

3.11 

99 

6986. xlO'3 

2.05 

3.136xl03 

4.01 

47 

1.149xl0'3 

3.98 

5.159xl04 

0.12 

4.  CONCLUSIONS 

In  this  paper,  we  have  adapted  the  structural  enrichment  meth  to  incorporate  non-linear  material 
properties  within  the  composite.  Using  structural  enrichment  requires  that  the  material  and  enriched 
shape  functions  be  updated  for  all  elements  at  each  time  step.  The  results  show  an  improvement  in  the 
thermal  strain  calculations,  for  the  example  shown,  when  using  the  structural  enrichment  techniques  over 
pure  homogenization,  ranging  from  0.5%  to  12%. 
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